#include "mloop.h"

#define INTERVAL 50.0
#define ROUND_MAX 150000
#define CHECKPOINT 5000

#define OUTPUT_THRESHOLD 0.05

int check_flat (int *H, float *S, float E_low, float E_high)
{
	int i, bin_low, bin_high, n, m;

	bin_low = E_low / INTERVAL;
	bin_high = E_high / INTERVAL;

	n = 0;
	for (i = bin_low; i <= bin_high; i++)
	{
		n += H[i];
	}

	m = 0;
	for (i = bin_low; i <= bin_high; i++)
	{
		if (H[i] * (bin_high - bin_low + 1) < n * 0.8)
			m ++;
	}

	if (m * 20 > bin_high - bin_low + 1)
		return 0;
	else
		return 1;
}

/* return : 1 converged */
int mmc(int max_output, float threshold)
{
	int i, n, m, num, num1, flag, flag1, flag2;
	int niter, norm, bin0, bin1, *H;
	double random_number;
	float *S, f, E_max, E_low, E_high, E_cut, E_cut2, converged, threshold0, threshold1, threshold_step, rmsd_avg;
	static float Converged = 9.9;
#ifdef DEBUG
	float rmsd;
	static float rmsd_min1 = 9.9, rmsd_min2 = 9.9, rmsd_min3 = 9.9;
#endif

	printf ("Max Output: %d\n", max_output);

	E_cut2 = loop.resnum * 2000.0 * OUTPUT_THRESHOLD;

	if (loop_length > 10)
		threshold_step = 0.3;
	else if (loop_length > 6)
		threshold_step = 0.25;
	else
		threshold_step = 0.20;
	threshold1 = threshold * 0.5;
	converged = 9.9;

	norm = 1000;
	H = (int *) calloc (norm, sizeof (int));
	S = (float *) calloc (norm, sizeof (float));
	E_max = norm * INTERVAL;
	E_cut = loop.resnum * 1000.0 * OUTPUT_THRESHOLD;
	for (i = 0; i < norm; i++)
		S[i] = 0.0;

	loop.circle = (2 * (loop.resnum - 1) > 10) ? 10 : 2 * (loop.resnum - 1);
	E_low = E_high = current[0].energy;

	flag1 = flag2 = 0;

	for (niter = 1; niter <= 100; niter ++)
	{
		f = 1.0 / niter; num = num1 = 0;
		for (i = 0; i < norm; i++)
			H[i] = 0;

		rmsd_avg = 0.0; m = 0;
		for (n = 0; n < ROUND_MAX; n++)
		{
			memcpy (&current[1], &current[0], sizeof (struct candidate));

			random_number = drand48 ();
			if (random_number < loop.probility)
				purturb_mainchain ();
			else
				purturb_sidechain ();

			if (current[1].energy > E_high)
			{
				E_high = current[1].energy;
				if (E_high > 3000.0 * loop.resnum)
					E_cut = loop.resnum * 3000.0 * OUTPUT_THRESHOLD;
				else if (1000.0 * loop.resnum < E_high)
					E_cut = E_high * OUTPUT_THRESHOLD;
			}
			else if (current[1].energy < E_low)
				E_low = current[1].energy;

#ifdef DEBUG
			rmsd = calc_rmsd (&ref, &current[1]);
			if (rmsd < rmsd_min1)
			{
				rmsd_min1 = rmsd;
				printf ("E1 %f\tE2 %f\tRMSD %f\n", current[1].energy1, current[1].energy2, rmsd);
			}
#endif

			if (niter > 1 || n >= 50000)
				flag = 1;
			else
				flag = 0;
			if (current[1].packing <= 2 && ((flag == 0 && current[1].energy2 < E_cut2) || (flag == 1 && current[1].energy2 < E_low + E_cut)))
			//if (current[1].packing <= 2 && (current[1].energy2 < E_cut2 || ((niter > 1 || n >= 50000) && current[1].energy2 < E_low + E_cut)))
			{
#ifdef DEBUG
				if (rmsd < rmsd_min2)
					rmsd_min2 = rmsd;
#endif
				if (current[1].packing == 0)
					threshold0 = threshold;
				else
					threshold0 = threshold + threshold_step;

				memcpy (&candidate[nsaves], &current[1], sizeof (struct candidate));
				flag = loop_cluster (threshold0);
				if (flag > 1)
				{
					num ++; numofsamples ++;
				}

#ifdef DEBUG
				if (flag == 2)
				{
					if (rmsd < rmsd_min3)
						rmsd_min3 = rmsd;
				}
#endif
				if (flag > 1 && num % CLUSTER_THRESHOLD == 0)
				{
					printf ("Samples: %d Nsaves: %d\n", num, nsaves);
					if (nsaves == nsaves2)
					{
						reinit (&current[0]);
						flag1 ++;
						if (flag1 == 3)
						{
							converged = threshold;
							if (converged <= threshold1 || (converged <= Converged && nsaves >= max_output * 0.75))
							{
								flag2 = 1;
								break;
							}
							else if (threshold >= 0.1)
							{
								threshold -= 0.1;
								printf ("Threshold1 decreased to %f\n", threshold);
							}
							flag1 = 0;
						}
					}
					else
					{
						flag1 = 0;
					}
					nsaves2 = nsaves;
				}
				if (nsaves == max_output)
				{
					flag2 = 1;
					break;
				}
			}

			if (current[1].energy > E_max)
			{
				norm += 100;
				H = (int *) realloc ((void *) H, norm * sizeof (int));
				S = (float *) realloc ((void *) S, norm * sizeof (float));
				for (i = norm - 100; i < norm; i++)
				{
					H[i] = 0;
					S[i] = 0.0;
				}
				E_max = norm * INTERVAL;
			}

			bin0 = current[0].energy / INTERVAL;
			bin1 = current[1].energy / INTERVAL;
			if (S[bin1] < S[bin0] || drand48 () < exp (S[bin0] - S[bin1]))
			{
				H[bin1] ++;
				S[bin1] += f;
				memcpy (&current[0], &current[1], sizeof (struct candidate));
			}
			else
			{
				H[bin0] ++;
				S[bin0] += f;
			}

			if (n >= CHECKPOINT && n % CHECKPOINT == 0)
			{
				if (check_flat (H, S, E_low, E_high))
					break;

				if ((niter > 1 || n >= 50000) && num == num1)
					reinit (&current[0]);
				num1 = num;
			}
		}
		printf ("Iterate: %d E_low: %f E_high: %f E_cut: %f Sampled: %d Num: %d Nsaves: %d\n", niter, E_low, E_high, E_cut, n, num, nsaves);

		if (flag2) break;
	}

	if (converged < 5.0)
	{
#ifdef DEBUG
		printf ("Length is %2d Converge to %.2f Best RMSD %8.3f %8.3f %8.3f\n", loop_length, converged, rmsd_min1, rmsd_min2, rmsd_min3);
		printf ("Total %d %ld conformations are saved\n", nsaves, numofsamples);
#else
		printf ("Converge to %.2f\tTotal %d %ld conformations are saved\n", converged, nsaves, numofsamples);
#endif
		if (converged < Converged)
			Converged = converged;

		flag = 1;
	}
	else
	{
#ifdef DEBUG
		printf ("Length is %2d Unconverged Best RMSD %8.3f %8.3f %8.3f\n", loop_length, rmsd_min1, rmsd_min2, rmsd_min3);
		printf ("Total %d %ld conformations are saved\n", nsaves, numofsamples);
#else
		printf ("Unconverged.\tTotal %d %ld conformations are saved\n", nsaves, numofsamples);
#endif
		flag = 0;
	}

	free (H);
	free (S);

	return flag;
}
